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Abstract. 
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exact associated Nakajima-Zwanzig master equation for quantum open system 
evolution. A mean field approximation for the memory kernels is introduced 
that yields, for an optimally chosen projection operator, a completely determined 
inhomogeneous master equation. Previous proofs of positivity and equilibration 
arc extended to these new inhomogeneous non-Markovian master equations. We 
study a nitrogen vacancy center in diamond interacting with 13 C impurities to 
illustrate the method. 
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1. Introduction 

Non-Markovian generalizations of the well known Kossakowski-Lindblad master 
equation^ have been a focus of many recent studies[2 l3ll4ll5ll6]l7ll8ll9lH0 l lTT l 
[T2l [TBI HH [TSl EH H7J [IB] in quantum open systems theory [JJj]. Such structures can 
preserve complete positivity[2j [3j [4[ [5l [6l [7] or positivity HUTU]. Efforts to quantify 
the non-Markovian character of such dynamics [16j [TTl [18] are also underway, as are 
explorations of relations between integrodifferential and time local forms [TTJ [T3]. The 
interest stems partly from the fact that equilibration in a Markovian model cannot 
simultaneously be complete positivity preserving and correct in the long time limit 
except in certain exceptional cases [5UJ[3T]. However, it is also true that non-Markovian 
effects are expected to play an important role in the low temperature condensed phase 
architectures [19] proposed for things like quantum computers [22j. 

Under certain conditions positivity preserving non-Markovian generalizations 
of the Kossakowski-Lindblad master equationpQ may also be shown to correctly 
equilibrate [8 . As with the Kossakowski-Lindblad master equationfTJ it is important 
to connect the parameters and memory functions of these abstract mathematical 
structures to data for specific physical systems. Semi-classical methods [23] and mean 
field approaches[M] have been explored. Here we examine a modification of the mean 
field method to connect a specific positivity preserving master equation which correctly 
equilibrates [5] to data for an electronic spin 5=1 associated with a Nitrogen- Vacancy 
center in diamond interacting with the 7=1/2 nuclear spins of 13 C impurities. 

We start in Section 2 with a form of the exact Nakajima-Zwanzig master 
equation[25] (NZME) which is valid for cases in which the bath operator A, associated 
with the NZME projection operator 

P=l s ® Atr B {-}, (1) 

commutes with the bath Hamiltonian. [Is is the identity operator in the system 
Hilbert space. Note that we will use Tr^ to denote a full trace over bath degrees of 
freedom and tr^ to denote the partial trace.] 

An ensemble of A operators is defined in Section 3 wherein each A is specified by 
setting the values of a set of scalar parameters. The NZME is exact for every parameter 
set. We then define a mean field approximation for this ensemble. Because the mean 
field master equations are approximate each will have a different accuracy. Most of 
the mean field master equations also possess an inhomogeneous term not previously 
considered in the abstract theory [5]. With the exception of Ref. [TU] previous 
studies of non-Markovian structures [71 151 l9l [TO] did not consider inhomogeneous cases. 
Unfortunately the sufficient conditions discussed in Ref. [10] are inapplicable here. 

In Section 4 we thus reexamine the issues of positivity and equilibration for 
this more general inhomogeneous structure. We find that it is sufficient for the 
memory functions to satisfy a set of constraints for positivity and equilibration to 
be guaranteed. Section 5 discusses the necessary basic structure of a model memory 
function that is consistent with the dictates of Section 4. In Section 6 we select the 
scalar parameters of the projection operator and those of the abstract theory such that 
the memory functions of the two different theories are optimally matched. Through 
this synthesis we are able to determine specific values of all parameters for the structure 
which preserves positivity and correctly equilibrates. 

In Section 7 we discuss the results obtained by applying the synthetic theory 
to a simplified model of an NV center in diamond. The predictions of the synthetic 
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master equation (SME) are compared to exact calculations obtained by well tested 
techniques 26 . Good accuracy is obtained at short and intermediate times but 
problems are seen in the long time limit. The diagonal elements are accurate within 5 
% but the off-diagonal elements, which in the exact case have persistent oscillations, 
are inaccurate. We discuss possible modifications of the abstract structures to correct 
this problem in Section 8. 



2. Nakajima-Zwanzig master equation 

Consider a simple composite system-bath with total Hamiltonian 

H tot = H®l B +S®B + ls®H B . (2) 

where H and Hb are the system and bath Hamiltonians, S and B are system and 
bath coupling operators, and H and S are in general non-commuting. Here Tb is the 
bath identity operator. We will employ units such that energies are measured in GHz 
and time is measured in ns. S will be unit-less and H, B and H B will have energy 
units. 

We define a projection operator P via ([1]) and choose the A operator for the bath 
with the following properties: 

At = A (3) 

Tr B {A} = 1 (4) 

[H B ,A] -0. (5) 

It is important to note that A ^ ps in general where ps = e~ HB l kBT /Zb is the 
canonical density operator and Zb = Trs{e~ ff£! / fcBT }. We also define a second 
projection operator via 

Q=I B -P (6) 

so that P + Q=1 B , P 2 = P and Q 2 = Q, and PQ = QP = 0. 

The total density operator ptot{t) — e~ iHtott ptot{0)e lHtott obeys the quantum 
Liouville equation 

dptot{t) 



dt 

where L to t = [Htot, •]• It can then be shown that 

dPptotjt) 
dt 

dQptotjt) 
dt 

and solving Eq. © for 



iL to tptot(t) (7) 
;n be shown that 

iPL tot Pp tot (t) - iPL tot Q Ptot (t) (8) 



= - iQL tot Pptot(t) ~ iQL to tQptot(t) (9) 



Qptot(t) = e-^ ito **Qp tot (0)-i / dt' e-^^-^QL^PpUt'^W) 

Jo 

and substituting back into Eq. ([8]) gives 

dPptot{t) = - iPLtote-^^QpUO) - iPL tot P Ptot {t) 
dt 

- [ dt' PLtotQe-^-^-^QLtotPptotit'). (11) 
Jo 



We choose an initial state of the form 

PtotiP) =p(0)®p B (12) 
with pure initial system state 

p(0) = |V(0))(V(0)|. (13) 



Note that 



where 



Pptot(t)= p(t)®A (14) 



p(t)=tr B { ftot (t)} (15) 
is the usual reduced density matrix of the system. In consequence 

PLtote-^^QptotiO) = PW <OW Qp(0) ® p B (16) 

= [5, ti B {Be- iQLMt p(0) ® (p B - A)}] ® A (17) 

PL tot Pptot(t) = PL tot p(t) <giA=[H + SB, p(t)} ® A (18) 

QLtotPptot (t) = QL t otp{t) ® A = [5(5 - B) , p(t) <g> A] , (19) 

where £? = Ttb{-BA}, and it follows that obeys an ensemble of master equations 
each of the form 

*Ejj!± = -i[H + SB,p(t)} i[S, (M^(t) - M^(t))p(0)} 

dt' [S, {^{M^\t-t') + M {i \t-t'))-BM^{t-t') \ [S,p{t')\ 

+ i (m< 3 >(< - f ) - M< 4 )(( - i')) [5, p(t')]+] (20) 
where [S,p(t')]+ = Sp(t') + p(t')S and 

M {1 \t) = tr B {Be- iQLtatt p B } (21) 

M {2 \t) = tiBlBe-^^A} (22) 

M< 3 >(t) = tv B {Be- iQLioit BK} (23) 

M<*>(t) = tr B {Be- iQLtott AB} (24) 

are memory operators. 

It is doubtful that any choice of A would lead to a set of (t) which have simple 
analytic forms. Numerical evaluation of M^ k \t) would be as difficult as solving the 
full system plus bath Schrodinger equation. Hence Eq. (|2H1) is exact but essentially 
unsolveable. 

3. Mean field approximation 

A solvable approximate master equation can be obtained by replacing the system 
operators M^ k \t) by their scalar mean values 

N 2 

< M(fe) W) = X>s{x}M<*>(*)^} (25) 

3=1 

for k = 1, . . . , 4 where Xj are a complete set of states for the system Liouville Hilbert 
space and N is the dimension of the usual state space. Obviously this is an uncontrolled 
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approximation that introduces an unknown amount of error. Nor is the approximation 
unique as we will show in Section 8. 

In this mean field approximation Eq. (I20[) yields solvable master equations 

^ = -i[H + SB, p(t)] - ^ (0) (t) [S, p(0)] 

- f dt' [S,K^(t-t')[S )P (t')} + K^(t-t')[S,p(t')} + } (26) 
Jo 

where 

Jf(°) (t) = (M ^ (t)) - (AfW (t)) (27) 
JjfW(t) = I((Af( 3 )(t)) + (AfW(t)» -F <AfW(t)) (28) 

X( 2 )(t) = ^«M( 3 )(t)>-(M( 4 )(t))) (29) 

are scalar memory functions that differ depending on the choice of projection operator 
A. Hence, the error will also be dependent on A. 

3.1. Self- consistent Born approximation 

We will employ an operator A, associated with the projection operator (JTJ, of the form 



A = P b+J2^( H b- H b)pb, (30) 



which clearly satisfies conditions j3j-([5]). The mean values of the operators M^ k \t) 
can then be approximated [24 using techniques from Random Matrix Theory and a 
second order self-consistent Born approximation [2T1 128] . The details are discussed in 
Ref. [23]. One subtlety is that the W(t) memory function of Ref. [24] has the property 
W(0) = 1 but this is not true of our (M^(t)). To facilitate application of techniques 
from Ref. [24] we define an operation such that each of Eqs. (|2i p -(|24 p can be written in 

the mean as (Af( fc )(£)) = e~ l Q Ltott k - Application of this identity for the total identity 
operator X = Xs ®Xb gives for k = 4, for example, (tiB{BXs ® XbAB}) = X 4 . Our 
memory functions can be handled according to the methods of Ref. [21] provided we 
define W(z) = (M^ (z)} /X k , as the Laplace tranform of W(t), where each k is treated 
separately and where (M^(z)) is the Laplace transform of (M^ k \t)). 
The application of methods from Ref. [24] then yields 

(MM(t))=BW(a 1 ,0 1 ,t) (31) 

(M^(t)) =BW(a 2 ,p 2 ,tJ_ (32) 

(M^(t)) = (MW(()) =W W(a 3: (3 3 ,t) (33) 
where B = Tr B {Bp B }, B 1 = Tr B {B 2 A}, and 

%,&,t) = T V (-a k tr( — ) J n/2+1 (fat) (34) 



where 



a k = (AAh - AA k )/yX k AAh (35) 



h = (AAh + AAk)Nl k AAh- (36) 



k 


2fc 


AAh 


AA k 


1 

2 
3 


B 
B 
B 2 


(tr B {BAA^p B }) 
(tr B {BAA*A}) 
(tT B {BAA*BA}) 


(tr B {BAA PB }) 
(tT B {BAAA}) 
(tr B {BAABA}) 



Table 1. Parameter formulas for the mean field approximation 



[Note that Ref. [53] employs notation like (A/4J) for an average over system and bath 

degrees of freedom instead of our notation *4„4t fe which in general does not denote a 
full average. Only the system part is a true average in our case.] In Eqs. (|35|) and (|36| 

A = QL tot and the formulas for these parameters AAk, AA^k aud Ik for different 
k are given in Table [T] Note again that Ik is a normalization constant introduced 
because (M^(0)) ^ 1. 

Note also that Eq. (|33|) implies that 

K^(t)=0. (37) 

Consequently the mean field master equation reduces to the simpler form 

^ = -t[H + SB, p(t)} - iKW (t) [S, p(0)] 

dt' K^(t-t')[S,[S,p(t')]}. (38) 



For k = 1 one can show that the general definitions of Table 1 yield the explicit 
forms 

TPAAi = {[H, [#, •]]) (S - B) + ([5, [JJJ) _0f - F - B(B - B)) 
+ ([H, [S, ■]]) {B 2 - SB) + ([S, [S, •]]) (5 3 - B 2 B - B£ 2 + £ 2 £)(39) 

h 2 AA\ = ([H, [H, ■]]) (B-B-A (Tr B {B} - Tr B {Z B }£)) 

+ (([S, [H, •]]) J- ([ffJ5, ■]]) )(W-BB- A(Tr B {5 2 } - Tr B {B}B)) 

+ ([£, [5, •]]) (B3_- 5 2 S - A(Tr B {5 3 } - Tr s {i? 2 }£)) (40) 

where caligraphic averages like £? 3 are computed using A, i.e. B 3 = Tr B {£> 3 A}, while 
others like B 2 are computed using ps, i.e. B 2 — Ttb{B 2 pb}- For k = 2 

ft 2 X4 2 = {[H, [S, ■]]} (F-fi 2 ) + ([5,[^-]]) (W-2Bm+B 3 ) (41) 

ft 2 A4t 2 = ([if, [H, ■]]) Tr^A 2 }^^}/? - Tr B {B}) 

+ (([S, [H, ■}}) _ + (mS, •]]) ){W -B 2 - Tr B {A 2 }(Tr B {f? 2 } - Tr B {B}B)) 
+ ([S, [S, •]]) (B 3 - B 2 B - Tr B {A 2 }(Tr B {S 3 } - Tr B {B 2 }B)) (42) 



and for k = 3 



ft 2 A4 3 = ([/f, [if, •]]) (£ 2 - B 2 ) + ([S, [if, •]]) (B 3 - 2B 2 B + B 3 ) 
+ ([H, [S, •]]) (B 3 - BW) + ([S, [S, •]]) (W-Wb-W 2 + b 2 W) 

+ Tr B {BH 2 B BA} - 2Tr B {H B BH B BA} + Ty b {H 2 b B 2 A} (43) 

h 2 AAh = ([H, [if, •]]) (g*- B 2 - Tr B {BA 2 }(Tr B {B} - Tt b {Ib}B)) 

+ ({[S, [if,-]]) _+([/^[5,-]]) )(W-WB-Ti B {BA 2 }(Tr B {B 2 }-Tr B {B}B)) 

+ ([S, [S, ■}}) (B 4 - B 3 B - Tr B {BA 2 }(Tv B {B 3 } - Tr B {B 2 }B)) 

+ Ty b {BH b BA} - 2Tt b {H b BH b BA} + Tr B {H 2 B B 2 A}. (44) 
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Note additionally that with projection operator (|30l) 

Hp 

W = ~W + Vj {B l H J B - ~B l H B ) (45) 
i=i 

n p 

A = PB + J2 IS (p bH b ~ PbH> b ) (46) 
j'=i 

and other quantities X can be similarly computed via X — Tt{Xpb} (e.g. ~pb = 
Tr B {p|}). 

All parameters for the mean field master equation are now well defined for a 
specific choice of A. We could in fact compare the results from master equation ([38]) 
to exact results computed using the methods of Reference [26]. The bath traces 
and averages in the above formulas, as in Ref. [26] . would be computed in the 
eigenbasis of Hb and would include only the first ub eigenstates. Thus, for example, 
Ttb{Ib} — tib- The parameter ng is determined by the temperature and it is assumed 
that states higher in energy are essentially unpopulated. 

However, master equation (|38|> may not preserve positivity and may not correctly 
equilibrate. Hence, we will instead try to match the time-scales of (|38p with a more 
abstract theory which does preserve positivity and docs equilibrate. We will primarily 
focus on matching the memory functions K^ l \t) for the two theories. The abstract 
theory is discussed next. 



4. Positivity requirements and equilibration 

We consider a master equation with the same general structure as Eq. (|3"8"]l but 
with potentially different memory functions. Our first goal, in Section 4.1, will be 
to show that memory functions exist for which the density matrix has a stochastic 
decomposition and therefore preserves positivity. By stochastic decomposition we 
mean that p(t) can be written as an average 

p(t)=M[Mt))m)\] (47) 
where M[-] is a mean evaluated over different stochastically evolving wavefunctions 
v(l))- 

Below in Section 4.2 we will examine the circumstances under which equilibration 
is possible. 



4--1- Stochastic decomposition 

Here we generalize the proof from Ref. [8] to the case of inhomogeneous master 
equations. We also relax one unnecessary restriction imposed in Ref. [5]. 

Following Ref. [S] we assume that we may write the Laplace transform K^(z) 
of the memory function K ^(i) in the form 

-mr^Km^^-^ < 48 > 

and we will require that V(t), the inverse Laplace transform of V(z), obey the 
inequalities 

V(t) > (49) 

V'{t) = dV{t)/dt < 0. (50) 



We will also assume that the effects of the drift governed by L = [H + BS, ■] and 
the dissipation governed by C — [S, [S, ■]} can be treated separately with the overall 
dynamics obtained via a Trotter product formula|10j. Thus, we will focus on the 
evolution equation 

dp(t) 



dt 



= -iK^(t)[S,p(0)\ 



dt' K {1 \t - t')[S 2 p{t') + P (t')S 2 - 2Sp(t')S] 



which after a Laplace transformation can be written as 

(zp(z) - p(0j) = -i ( R(z) + ^M) 



[S,P(0)] 



- [S 2 p(z)+p{z)S 2 -2Sp{z)S] 
where we have defined the Laplace transform R(z) of a function R(t) via 

Using Eq. (j48|) and the fact that (f5Tj) implies 



dp(t) 
dt 



iK(°\0)[S,p(0)}, 



wc then have 



z(zp(z) - p(0)) - -p(t) \ t=Q -X(zp(z) - p(0)) 



= -iK^(0)R(z)[S,p(0)] 

- V{z)(zp{z) - p(0)) - K^(0)[S 2 p(z) + p{z)S 2 - 2Sp(z)S]. 
It follows that we may write 



d 2 p(t) dp(t) 



dt 2 



dt 



l dt' V (t-t>) d J*n 

' dt 1 



(51) 



(52) 
(53) 

(54) 



(55) 



- iK^(0)R(t)[S,p(0)\ - K^(0)[S 2 p(t) + p(t)S 2 - 2Sp(t)S] (56) 

= -V(0)p(t) + V(t)p(0) - f dt' V'(t t')p(t') 

Jo 

- iK^{0)R{t)[S,p{0)} - K^(0)[S 2 p{t) + P {t)S 2 - 2Sp(t)S], (57) 
by inverse Laplace transformation of Eq. (|5"5j) . Noting that 



dt 1 at 1 dt 



(58) 



and substituting into Eq. (|57p we finally have 



d 2 [p(t) 



,-\t/2] 



dt' 2 



= [A 2 /4- V{0)]p(t)e- M/2 + V{t)e- xt/2 p{0) 

dt'V'(t - t')e- A ^- t,)/2 p(t')e- xt '' 2 - iK^^R^e-^iS, p(0)} 
K (1 \Q)[S 2 p{t)e- xt / 2 + p{t)e- xt ' 2 S 2 - 2Sp{t)e- xt l 2 S] (59) 



which has a stochastic decomposition provided that Eqs. P5|) - (l5U)) . A 2 /4 > V(0), as 
well as other constraints discussed below are satisfied. Note that Ref. [5] required 
that A 2 /4 = V(0) which we will see is unnecessary. 

To avoid dealing with specific initial conditions p(0) we employ the Hadamard 
representation[9]. The Hadamard representation of the propagator is defined via 
(s\p(t)\s') — U s , s '(t)(s\p(0)\s') where S\s) = s\s). If the eigenvalues of U S:S i(t) arc 
positive then p{t) will be positive semidefinite for all positive semidefinite initial 
densities^. The propagator can also be calculated from U S)S >(t) — M[u s (t)u*,(t)] 
using Eq. (14"T1) where (s\ip(t)) = u s (t)(s\ip(0)) and u s (t) is again defined in the 
Hadamard sense. Rather than deal directly with u s (t) and U s ^ s '(t) we defined instead 
v s (t) = u s (t)e~ xt / 4 and V s , s .{t) = U s . s ,(t)e~ xt / 2 where V s , s .{t) = M[v s (t)v*,(t)]. 

Let a t = — R(t)/\R(t)\ then provided that \R(t)\e~ xt / 2 is monotonically 
decreasing, the ltd equations [T^] 

dv s (t) = ia t ^K^HQ) sv s (t) dt + y s (t) dw { t 1] 
+ ^\R(t)\e- xt /^K( 1 ){Q)dw[ 2) (60) 
dy s (t) = - ia t ^KW{fi) sy s {t) dt + AW)^ Xt/4 dw[ 3) 



+ / d4 4 V-V"(*-*')e" A(t ~ t ' )/ V0 dw[ 5) 



+ ^-±(\R(t)\e^/^KW(0)y w l 6) 

+ ^X 2 /4-V(0)v s (t) dw [ t ] (61) 

yield the equations 

J t V S {t)v*At) = iaty/KWm* ~ S')v s (t)v* s ,(t) 

+ y s {t)y* sl {t) + \R{t)\e- xt / 2 ^K^(0) (62) 
j t Vs(t)y*At)= - ivty/KW(0)( s - S r )y s (t)y:,(t) + V{t)e- Xt ' 2 
+ f dt'[-V'(t - t'^e-^-'^^Vsit'X^t') 



+ (X 2 /4-V(0))v s (t)v* s ,(t) 

- | {\R{t)\e-^yjK^j (63) 

in the mean. Note that dw t are real independent Wiener processes [19]. Provided 
that the sign of R(t) does not change, taking the derivative of (|6"2"|) and using (1B"3"|) 
then gives the equation 

= {\ 2 A-V{0))V^(t)-^R{t)K (1 \^- M,2 {^^) 

+ V{t)e- Xt ' 2 - [ dt'V'(t - i')e _A(t_ *' )/2 V s ,^(t') 
Jo 

-KW(0)(s-s') 2 V s , s ,(t) (64) 
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which is the Hadamard representation of Eq. (|59")1 . Initial conditions v s (Q) = 1 
guarantee that V s s /(0) = 1. Hence, we find that 
d- 
dt 



jv s {t)v:,{t)\t=0 = ia^K^(Q)(s - s>) + y s (0)y* (0) 



+ \R(0)\^KW(0) (65) 
which we require to equal 

j t v 8 (t)v*,{t)\ t=Q = -iK^(0)(s - s') - A/2 (66) 

in order to match the original initial condition (|54[) . Thus, it appears that we must 
choose 

KW(0) = -a oy /KW(p) (67) 

-A/2 > |i?(0)|y / if( 1 )(0) (68) 

so that we can pick initial values y s (0) = J-X/2 - \R(0)\ sjK ^(Oj. This clearly 
implies that 

A < 0. (69) 

Since the Hadamard propagator has a stochastic decomposition it is positive 
semidefinite and so p{t) will be positive semidefinite. This completes the proof of 
positivity for pit) provided that a memory function satisfying (j^5|) - ([5Uj) and (1671) - (f6"8|) . 
A 2 /4 > 1^(0), and a t constant and \R{t)\e~ xt / 2 monotonically decreasing, can be 
found. We will discuss an appropriate set of memory functions in Section 5. 

4- 2. Equilibration 

The following argument generalizes results from Ref. [S] to the case of an 
inhomogeneous master equation. 

Consider a non-Markovian Kossakowsi-Lindblad^l (NMKL) master equation of 
integrodifferential form[7J [8l [9J [10] which employs a scalar memory function K(t) and 
has a single Kossakowski-Lindblad dissipation operator constructed from the system 
interaction operator S: 

^-=F(t)-i[H + BS,p{t)] 

dt' K(t- t')[S 2 p{t') + p(t')S 2 - 2Sp(t')S] . (70) 



Here F(t) represents the inhomogeneous term. 

For additional notational simplicity we will introduce operators L = [H + BS, ■] 
and C — [S, [S, ■]]. Laplace transforming equation (ffD)) yields 

zp{z) - p(0) = F(z) - iLp(z) - K(z)£p(z) (71) 
from which one can then show that 

p(z) = {z + iL + K{z)C)- 1 [p{Q) + F(z)} (72) 
where p(z) = dt e~ zt p(t) is the Laplace transformed reduced density matrix, 
K(z) is the transform of the memory function and F{z) is the transform of the 
inhomogeneous term. We then define an operator 

G{z) = {z + iL + k{z)Cy 1 (73) 
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so that p{z) = G(z)[p(0) + F(z)]. 

Defining an operator Gq(z) = (z + iL)^ 1 , and letting \n) and E n denote the 
complete set of eigenvectors and eigenvalues of H + BS, it then follows that 

Km zG (z)\n)(n\ = \n)(n\ (74) 

2 — >0 

lim zG (z)\n)(m\ = 0. (75) 

z— >0 

[We have assumed that the spectrum of H + BS is discrete. Generalization to the case 
where all or part is continuous is straightforward.] In fact lim 2 ^o zGo(z) — Ho where 
IIo = l n )( n |Trs{|n)(7i|-} is a projection operator. 
It then follows that 

G(z) = {G (z)- l + K(z)£)- 1 (76) 

= [G Q (z)-\l + G {z)k(z)C)]- 1 (77) 

= (l + G (z)K(z)Cy 1 G (z) (78) 
Hence from Eq. (|72|) we obtain 

p(oo) = lim zG(z)\p(0) + F(z)} (79) 

2— >0 

= lim(l + G (z)K(z)£)- 1 zGo(z)[ j o(0) + J F(z)] (80) 
= hm(l + G (z)K(z)C)- 1 U o \p(0) + F(z)] (81) 

-1 



lim (l + zG (z)(A-(z)/z)£) n o [p(0) + F(z)}. (82) 



Thus, if the limit 



exists and 



then 



li m = K (83) 

2^0 z 



lim F(z) = (84) 

Z^rO 



p(oo) = (1 + /dlo^IIoKO) (85) 

and p(po) satisfies IIop(co) = p(oo). Perhaps this is more clearly seen by trivially 
rewriting ([55]) as p(oo) = n o p(0) — Kil £(l + Kn o £)~ 1 IIop(0) and recalling that 
IIq = Hq. In addition, since IIop(O) = J2 n ( n \p(®)\ n ) l n )( n l * ne l° n § time limit depends 
at most on the diagonal elements of the initial density matrix. Condition (|84|) is 
necessary to prevent dependence on the off-diagonal elements of the initial density 
matrix. Thus partial equilibration is possible in this inhomogeneous NMKL model. 
Exact calculations for some model systems equilibrate in a similar manner [2^1 130j . 



5. SME memory functions 

To find appropriate memory functions K^°\t) and K^(t) we start with a fairly general 
polynomial model for K^ 1 ' (z) which obeys the basic constraints discussed above and 
in Ref. |S]. Specifically, we assume that 

KW{z) = KW{0) 3 I ■ W 

z a + pz z + vz + 7 
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It is necessary that the polynomial in the denominator be of one degree higher than 
that in the numerator so that a non-zero initial value K^(0) is guaranteed for the 
memory function. The factor of z in the numerator guarantees a non-zero limit k in 
Eq. (|83p provided that (3 and 7 are non-zero. This also guarantees that the memory 
function is integrable and so has a proper Markovian limit. For non-zero j3 either v 
or 7 must be non-zero to guarantee that memory is of finite duration. 
Given the memory function (|86p it then follows that we must have 



A = /3 — fx (87) 

V{z) = [v-^ + p 2 + 1 -]-±— . (88) 

z z + p 

The k of Eq. (|83p can also be obtained as k — (0)^/7. Now V{t) will be positive 
if V(z) is completely monotone [31]. If v — [i(3 + f3 2 > and 7 > then the first factor 
in Eq. ([88]) will be completely monotone. If f3 > then the second factor will also be 
completely monotone. The product of completely monotone functions is completely 
monotone [3T], Hence, V(z) is completely monotone and constraint (|49| is satisfied, 
i.e. V{t) > 0. 

Note also that 

WO) = lim zV{z) = v-^I3 + P 2 . (89) 

z— >oo 

Denote A(<) = — V'(t), which must be positive by ([50]) . and its Laplace transform will 
be A(z) = V(0) - zV(z). Inserting Eq. ([88]) gives 

z + p 

which is clearly completely monotone and so A(t) will be positive. 

For K(°\t) we know that (H7])-([nH]) must be obeyed (see section 4.1), R(t)e~ xt/2 
defined via ([53"]) must not change sign and should be decreasing in magnitude (see 
section 4.1), and additionally (l84l) must be satisfied. The model 

RW{z) = g(°)(0) , Z[ \ +01) (91) 



satisfies all of these constraints. Clearly it follows from definition ([53]) . as well as 
and ([Hi]), that 



#(°)(0) 



Q 



so that 



and 



= 1^ - V (92) 

fl(t )e-^ = ^j|i e-i^ (94) 

which is monotonically decreasing in magnitude provided that 

3/3 > fi, (95) 
and it does not change sign. 
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Note also that 

lim K<V (z) = lim ( ) z(z + q) (96) 

z^Q z^Q Z A + /iZ 2 + VZ + 7 

= 0, (97) 

and since F(z) = -i[S, p(0)]K (0) (z) it follows that flU} is satisfied. 
Finally, conditions (p7)) - (|55|) require that 



\K (0 \0)\ = y/KW(0) (98) 

P > a > /3 + A/2. (99) 

We arbitrarily set a = j3 + A/4 and K^°\0) = \J KWjjfy. In the calulations reported 
below K^°'(t) is insensitive to the value chosen for a, and the reduced density matrix 
elements are insensitive to K^°'(t), 

6. Nitrogen- Vacancy model 

The native Hamiltonian is that of an S = 1 electronic spin associated with a single NV 
impurity center in diamond in its ground electronic state, and it takes the form|32] 

H = h x S x + DS X + E(S 2 X - S Y ) (100) 

where the magnetic field is along the ^-direction, h x = g e /i e B x , and Sx, Sy and Sz 
are operator components of the electron spin with the form|32| 

1 
S x = -5= ( 1 1 
1 




S z = 

The bath consists of 18 13 C impurity nuclear spins with Hamiltonian [33] 

17 18 

ff B = /4°)]T^) + (/?/C)£ Yl C jlk (3I<p4Q -I®.*")) (101) 

and[3J] /ii 0) = gnHnB x , 

C jtk = [1 - 3( Zj - zk f/lr, - r fe | 2 ]/| rj - r fe | 3 , (102) 

where C 2 = 5Z^=i Ejt=j*+i ^j,k- Here Zx , Jy and denote the components of 
the vector spin 1^ for the jth nucleus, and have the form of their corresponding Pauli 
matrices multiplied by I = 1/2. The factor of 1/C is required to correct the units 
and guarantee a proper thermodynamic limit. The impurity locations were selected 
randomly from a spherical integer lattice with a radius of 5 units. This means about 
4 % of the lattice sites are occupied by 13 C impurities whereas in real diamond the 
natural abundance is 1.1 % [35]. Of course we should really choose the correct diamond 
lattice locations|32| but our model has more serious flaws. 
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NV model Parameter 


Value (in GHz) 


Reference 


K 


.194 


m 


D 


2.88 


[32] 


E 


.1 




A X x 


.2 


m 


A X Y 


.02 


m 


A X z 


.02 


m 




1.08 x 10" 3 


m 


P 


4.52 x 10~ 5 


m 


k B T 


.0003 


36 



Table 2. Numerical values of NV model parameters for B x = 11 G 



Mean field memory function parameter 


Value 


B 


9.3513 x 10" 2 


B 


9.3276 x 10~ 2 




9.8692 x 10^3 


ai 


1.4111 


01 


1.4259 


CC2 


1.3935 


02 


1.3951 


a 3 


1.1953 


fa 


1.7843 



Table 3. Numerical values of parameters for memory functions of mean field 
master equation 



The system-bath coupling should consist of a complete tensor coupling all system 
spin components with all bath spin components [32], but the present formulation of the 
master equation can only handle a single system coupling operator. We thus choose 
arbitrarily to include only the x-components, i.e., we pick a system-bath coupling 
operator of the form 

H SB =S X ® (l/A) A k (A X x4 k) + A X Y^ k) + A xz I {k) ) (103) 

k 

where [33] 

^ fe = (l-3z 2 /|r fe | 2 )/|r fc | 3 (104) 

and A 2 = Ek=i A l Hence, our model is not completely correct at present. Until 
we generalize the abstract theory to multiple system coupling operators this is the 
best we can do. The total Hamiltonian is thus of the form @ with S = Sx and 
B = {l/A) J2 k Ak(A X xIx + Axyly + A X zli k) ) with parameter values given in 
Table 2 for a magnetic field of 11 G[521|55). 

The initial density matrix was chosen to take the form (1131) with 

|^(0)> = i=[V3|->+ l |0> + |+>] (105) 

where |— ), |0) and |+) denote eigenstates of Sx with corresponding eigenvalues — 1, 
and 1. 
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SME parameter 


Value 




45.9675 




46.4375 


V 


21.6505 


7 


106.1616 


B 


9.3276 x 10~ 2 




1.1665 x 10" 3 



Table 4. Numerical values of parameters for memory functions of SME 



6.1. Determination of parameters of SME 

Equations (|39)) -(|44 | involve quantities 

([H, [H, •]]) = (D + E) 2 /3 + 4h 2 x /3 + (D - 3E) 2 /9 (106) 

([H, [S, ■]]) = 4h x /3 (107) 

([S, [H, ■]]) = Ah x /3 (108) 

([S, [S, ■}}) = 4/3 (109) 

which can be evaluated given S and H for the NV model. All other parameters 
including the rjj of (|3"0|) for j = 1, • • • , n p and the /3, /i, v and 7 of (|5d| were set during 
a simulated annealing 35^ process where the parameters were optimized so that the 
mean field memory function K^\t) of Section 3 and the formal mathematical model 
of Section 5 match as well as is possible. We picked n p = 10 and required 1^1 < 300. 
For the parameters of K^\t) we defined variables Xi ■ ■ ■ , X4 via 

= X! (110) 

fi = /3 + X 2 + 2^fx~ 3 (111) 

v = X z + ^P-p 2 (112) 

7 = X 4 (113) 
and we required < Xk < 200. The function 

f^, . . .,mo,*i, ■■■,X 4 )= dt' [K$ F (t') - K^ IE (t')} 2 , (114) 

Jo 

where K^ F {t) is the mean field memory function calculated via ([3"T j) - (|55|) . and 

K { s ] ME {t) is the SME memory function based the polynomial model of Section 5, 
was minimized using the program simann.f of Ref . [35] . 

The optimal memory function parameter values are given in Table 3 for the 
mean field master equation and Table 4 for the SME. Note that conditions 3/3 > /i, 
13 > 0, 7 > 0, V(0) = v - n/3 + (3 2 = 4.5775 x 10~ 2 > 0, A = -0.4700 < 0, and 
A 2 /4 - V(0) = 9.45 x 10~ 3 > are all satisfied. 

The resulting mean field K^ 1 ' (t) (solid curve) and the polynomial model (dashed 
curve) are compared in Figure 1. The agreement is quite good. Note that K^\t) 
decays to zero after about 20 time units or about 2 x 10~ 8 s. We will however be 
exploring dynamics for times as long as .1 ms. One might be tempted to conclude 
that a Markovian formulation would work for this problem. However, we know [20ll21j 
that the Markovian form cannot capture equilibration correctly in general unless the 
coupling operator takes a very restrictive form. 



16 




Figure 1. KW(t) vs. t 

Our model for K^°'(t) (dashed curve) shown in Figure 2 is roughly similar to 
its mean field counterpart (solid curve) except at very short time. A higher order 
polynomial model would probably improve the agreement. 

6.2. Numerics 

The exact equations were obtained using the methods of Ref . 26 . A temperature of 
k B T = .0003 GHz, or T = 1.44 x 10~ 5 K, was chosen for which a total of n B = 20 
bath eigenstates make significant contributions. (Temperatures an order of magnitude 
lower can be achieved in practice |3B].) All exact calculations and SME parameter 
calculations were carried out in this eigenbasis. 

The SME (|38|) can be mapped to a set of four differential equations by finding 
the roots of the polynomial z 3 + pz 2 + vz + 7 = 0. From the three roots Zk one can 
find K (1 \t) = i2l =1 a k e Zkt . Then defining tt k {t) = J*dt' a k e* k ^-^[S,[S,p{lf)]] we 
obtain 



dp(t) 



dt 
dt 



= -i[H + SB,p(t)]-iK 



(0) 



(t)[s,p(o)]-Y,sik(t) 



k=l 



a k [S, [S,p(t)]]+z k n k {t) 



(115) 



(116) 



for k = 1, . . . , 3 with initial conditions f2fc(0) = 0. All ordinary differential equations 
were solved using standard Runge-Kutta algorithms [551 13"5] . 

7. Results 



The master equation and exact results are not visibly different until quite late in 
the dynamics. Figure 3 shows the mean spin components S x (t) — Tr{Sxp{t)}, 
S y (t) — Tr{SV/f(i)}, and S z (t) = Tr{Szp(t)} computed exactly (solid curve) and 
using the master equation (dashed curve) on the time interval [1000, 1050] (our time 
units are ns) . By this point some discrepancy is discernible. The amplitudes of the 
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0.04 ,- 
0.03 ; 



0.02 -i 




-0.03 1 ' ' 1 1 1 1 

5 10 15 20 25 30 

t 



Figure 2. K^(t) vs. t 

oscillations in S x (t) are diminished compared to those in the exact calculations and 
they are slightly phase shifted. The results for S y (t) are still in good agreement with 
the exact calculations. In the case of S z (t) the error is larger but it appears mostly as 
a shift of the overall envelope of the oscillations. 

Figure 4 shows the mean spin components on the time interval [4000,4050]. By 
this time the master equation results for S x (t) and S y (t) are further out of phase with 
their exact counterparts, and both amplitudes are diminished. For S z (t) the master 
equation results are completely out of phase with the exact results but the amplitudes 
are not any worse than they were at the earlier time. 

Figure 5 shows the purity V(t) = Trs{p(t) 2 } plotted on a much longer time 
interval [0, 100000] (or to .1 ms). The solid curve is the exact result while the 
dashed is that of the SME. Here it is clear that there are serious problems with the 
predictions of the SME. A purity of about .8 should be seen at long times but the 
SME predicts just .65 which it out by nearly 20 %. 

To try to trace the origin of the errors we examined the diagonal elements of the 
reduced density operator in the eigenbasis of H + BS. These are shown in Figure 6 
(a) for state 1 which is the lowest in energy, and for state 3 in (b) which is the highest 
in energy. The solid curves denote exact results while dashed curves indicate SME 
predictions. The relative errors in the lowest energy matrix element pi,i(t) are on 
the order of .1 %. The errors in the highest energy matrix element P3,3(i) are on the 
order of 5 %. The middle matrix element (not shown) is exact and constant because 
this eigenstate is also an eigenstate of Sx- Hence, the diagonal elements are not the 
primary source of the errors seen in Figure 5. 

Figures 7-9 show the real and imaginary parts of <Ji^(t), 0"i,3(i) and 02,3^) 
respectively in the eigenbasis of H + BS where a(t) = e K H +BS)t p (^ e ~i(H+BS)t ig 
the reduced density matrix transformed to a rotating frame. Solid curves denote 
exact results while dashed curves are the SME predictions. Figure 7 for (J\^{t) shows 
decent agreement between the SME and exact results until around 10000 time units. 
However, the master equation results decay to zero by about 40000 time units, while 




Figure 3. Mean spin matrices on [1000, 1050] 
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Figure 5. V{t) vs. t 

the exact results continue to oscillate with large amplitude. Similar results are seen in 
Figure 8 for ai. 3 (t). Figure 9 for er 2j 3(i) shows terrible agreement between the SME 
and exact results at all times except t = 0. 

The primary discrepancy between the master equation predictions and exact 
results thus arises in the off-diagonal elements in the long time limit. Moreover, 
the corresponding discrepancy in the purity suggests that these errors arise from 
problems with the dissipative terms in the SME. Since the results change little when 
we arbitrarily set K^°\t) = we also know that the inhomogeneous term is not the 
source of error. 

7.1. Numerical error 

We noticed an unexpected error arising in the eigenvalues of the density matrix 
computed using the SME master equation. One eigenvalue should remain zero at 
all times but we discovered that it becomes slightly negative at short times but 
then returns to zero. This problem persists in double precision even for a requested 
absolute and relative tolerance of 10~ 16 but vanishes in quadruple precision with a 
requested absolute and relative tolerance of 10~ 17 . Progamming in quad precision is 
impractical however since this greatly slows the computations. It does not seem likely 
that numerical errors of this type are responsible for the discrepancies between the 
exact and master equation results. 

8. Discussion 

The synthetic approach introduced in this paper yields a well defined master equation 
with no free parameters which preserves positivity and correctly equilibrates to a long 
time limit which commutes with the shifted native system Hamiltonian and which 
depends at most on the diagonal elements of the initial density in the Hamiltonian 
eigenbasis. Numerical tests for a model system consisting of an NV center in diamond 
interacting with 13 C impurities show good accuracy at short and intermediate times. 
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0.7786 




0.7772 1 ' 1 1 ' 1 

20000 40000 60000 80000 100000 

t 

(a) pi,i(t) vs. t 

0.0228 | 1 1 1 1 1 




0.0214 ' ' ' ' ' 1 

20000 40000 60000 80000 100000 

t 

(b) P3,s(t) vs - * 
Figure 6. Diagonal elements of density matrix 

The diagonal elements of the density in the shifted native Hamiltonian eigenbasis 
are also in good agreement with exact results at long times. There is however a 
serious problem with the off-diagonal elements. In a frame rotating with the native 
Hamiltonian the exact results show persistent oscillations with a period of about 10 
fis, while the master equation results decay to zero quite rapidly. 

One possible solution to this problem would be to consider forms of the projection 
operator A which do not commute with H& or B in which case the mean field 
master equation would have a structure similar to Eq. ([26]) but wherein the 
K^{t) memory function is non-zero. Clearly this structure will deviate from the 
standard Kossakowski-Lindblad form but there may be some circumstances under 
which positivity is still preserved. It is also clear that this would change the off- 
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(a) Re cti,2 (t) vs. t 



40000 60000 
t 



(b) Im (7i,2 (t) vs. t 
Figure 7. o"i,2(i) element of density matrix in rotating frame 



diagonal elements but leave equilibration unaffected. 

Alternatively, one might consider a more complex structure to accompany the 
modified A. Consider for example that one of the kernels of Eq. (|20|) has the form 



i(M( 3 )(t)-MW(t))=itr B {S[e 



—iQLtatt 



,B]A} 



(117) 



which can be rewritten using Kubo's formula[3D] as 



dt' tr B {B, 



,-iQL tat (t-t') 



[QL tot ,B]e 



-iQLtott' 



A} 



(118) 



23 



0.15 



V v V V 



t 

(a) Re cti,3 (t) vs. t 



(b) Im (7i,3 (*) vs - t 



80000 100000 



Figure 8. o"i,3(t) element of density matrix in rotating frame 



which would suggest that K^ 2 \t) in Eq. (|26|) be replaced by some operator like 
K^\t)[H, •] + K^\t)[S, •] which would be of Kossakowski-Lindblad form if K^\t) = 
0. Here K^\t) would then probably need to be proportional to K^{t) in order to 
preserve positivity. 
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